Example: Quadratic Regression Model

For illustration, we will assume that the function \(\mu\) is a quadratic function of \(x\). So,

\[\mu(x) = \beta_0 + \beta_1 x + \beta_2 x^2\]

and

\[Y = \beta_0 + \beta_1 x + \beta_2 x^2 + \epsilon\]

The coefficients of the quadratic polynomial \(\beta_0, \beta_1, \beta_2\) are unknown parameters to be estimated. The parameter \(\beta_0\) is directly interpretable as the mean value of the response when \(x=0\). The other parameters \(\beta_1\) and \(\beta_2\) are not readily interpretable.

Somewhat confusingly, the quadratic regression model above is an example of a linear model. In a linear model, the unknown parameters, not the known predictors, must enter the formula linearly.

Assumptions

We make three important (and one less important) assumptions about the errors. First, we assume that \(E(\epsilon|x) = 0\). Second, we assume that the errors are independent. Third, we assume that the variance function \(\mbox{Var}(Y|x)\) is constant; the unknown constant variance \(\sigma^2\) is an additional parameter of our regression model. Errors are often assumed to be normally distributed, but, for most purposes, the normality assumption is not necessary.

To estimate the unknown parameters of the quadratic regression model, we observed \(y\) and \(x\) for \(N\) subjects. The dataset could be presented in tabular form as follows.

\[ \begin{matrix} y_1 & x_1 \\ y_2 & x_2 \\ \vdots & \vdots \\ y_N & x_N \\ \end{matrix} \]

We can write the quadratic regression model as \[y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \epsilon_i \hspace{3em} i = 1, 2, \ldots, N\]

Matrix Representation of the Linear Model

It will be more convenient to write the model in matrix/vector representation. The regression equation is then written as

\[ \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \\ \end{pmatrix} = \begin{pmatrix} 1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ \vdots & \vdots & \vdots \\ 1& x_N & x_N^2 \\ \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \end{pmatrix} + \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N \\ \end{pmatrix} = \begin{pmatrix} \beta_0 +\beta_1 x_1 + \beta_2 x_1^2 \\ \beta_0 +\beta_1 x_2 + \beta_2 x_2^2 \\ \vdots \\ \beta_0 +\beta_1 x_N + \beta_2 x_N^2 \\ \end{pmatrix} + \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N \\ \end{pmatrix} \]

or, more compactly, as

\[y = X \beta + \epsilon\] where

\[ y=\begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{pmatrix} \hspace{3em} X=\begin{pmatrix} 1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ \vdots & \vdots & \vdots \\ 1& x_N & x_N^2 \\ \end{pmatrix} \hspace{3em} \beta=\begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \end{pmatrix} \hspace{3em} \epsilon=\begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N \\ \end{pmatrix} \]

Example: Total Cholesterol as a Function of Age

To illustrate the use of linear models with a quantitative outcome variable, we will investigate the relationship between serum total cholesterol (mmol/L) and age (years) for the US adult (ages 20-79) population using the NHANES dataset. For illustration, we will use a quadratic regression model; in practice, we would substitute a natural cubic spline for the quadratic polynomial.

Many procedures in SAS can fit linear models for quantitative outcomes. We will use PROC GLMSELECT, which includes the EFFECT statement. Other commonly used procedures for linear models for quantitative outcomes are PROC REG, PROC GLM and PROC GENMOD.

LIBNAME NHANES "../../../BMI552_Datasets/NHANES_May2023/";
*LIBNAME NHANES "~/my_shared_file_links/u48718377/NHANES";
OPTIONS FMTSEARCH=(NHANES.formats_NHANES work) NODATE NONUMBER;
TITLE1; TITLE2;
ODS NOPROCTITLE;
ODS GRAPHICS / LOESSMAXOBS=20000;

DATA Analysis;
  SET NHANES.NHANES (RENAME=(Race1=Race));
RUN;

ODS EXCLUDE ALL;
PROC GLMSELECT DATA=Analysis;
  EFFECT AgePoly = POLYNOMIAL(Age / DEGREE=2);
  ODS OUTPUT ANOVA=A FitStatistics=F;
    
  MODEL TotChol = AgePoly / SELECTION=NONE;
    
  STORE QuadraticModel;
RUN;

PROC PLM RESTORE=QuadraticModel;
  ODS SELECT ParameterEstimates;
  SHOW Parameters;
RUN;
ODS EXCLUDE NONE;

The least squares estimates of the regression coefficients are \[\hat{\beta}_0 = 3.1126, \hat{\beta}_1 = 0.08122, \hat{\beta}_2 = -0.00075\] The estimated mean serum total cholesterol as a function of \(x\) (age) is \[\hat{E}(Y|x) = 3.1126 + 0.08122 x - 0.00075 x^2\] This estimate is only valid for ages within the range of the observed data. Applying this estimate outside of the range of the observed data is called extrapolation.

DATA Score;
  DO Age = 20 TO 80 BY 1;
    OUTPUT;
  END;
RUN;

ODS EXCLUDE ALL;
PROC PLM RESTORE=QuadraticModel;
  SCORE DATA=Score OUT=Score Predicted LCLM UCLM LCL UCL / ALPHA=0.03125;
RUN;
ODS EXCLUDE NONE;

DATA All;
  SET Analysis Score(IN=S);
    
  Score = S;
RUN;

PROC SGPLOT DATA=All NOAUTOLEGEND;
  SCATTER X=Age Y=TotChol / MARKERATTRS=(SYMBOL=CircleFilled COLOR=LightBlue);
*  BAND X=Age Lower=LCL Upper=UCL / FILLATTRS=(COLOR=Red TRANSPARENCY=0.8);    
*  BAND X=Age Lower=LCLM Upper=UCLM / FILLATTRS=(COLOR=Red TRANSPARENCY=0.6);
  SERIES X=Age Y=Predicted / LINEATTRS=(COLOR=Red);
RUN;

Estimating the Variance

Since the variance \(\sigma^2\) is essentially the average of \(\epsilon_i^2\), we should expect that its estimator could be obtained by averaging the squared residuals \(r_i^2 = \hat{\epsilon}_i^2\). The maximum likelihood estimator of \(\sigma^2\) can be obtained by dividing the residual sum of squares \(\mbox{RSS} = \sum_{i=1}^N{r_i^2}\) by the sample size \(N\). This estimate is rarely used.

An unbiased estimator of \(\sigma^2\), which is obtained by dividing the residual sum of squares \(\mbox{RSS} = \sum_{i=1}^N{r_i^2}\) by its degrees of freedom (df), where the degrees of freedom is the number of observations minus the number of parameters in the mean function. For our quadratic regression model, the residual df is \(N-3\), so the estimate of \(\sigma^2\) is

\[\hat{\sigma}^2 = \frac{RSS}{N-3} = \frac{\sum_{i=1}^N{r_i^2}}{N-3} = \frac{\sum_{i=1}^N{(y_i-\hat{y}_i})^2}{N-3}\]

This quantity is called the residual mean square. In general, a sum of squares divided by its degrees of freedom is called a mean square. The square root of the residual mean square \(\hat{\sigma}\), the root mean square is also called the standard error of the regression. It is an estimate of the standard deviation of the error distribution; it is in the same units as the response variable.

In our example, \(\mbox{RSS} = 7674\), \(\hat{\sigma}^2 = 1.061\) and \(\hat{\sigma} = 1.030\).

PROC PRINT DATA=A NOOBS LABEL;
  WHERE Source="Error";
  VAR DF SS MS;
  LABEL DF="Residual df" 
        SS="Sum of Squares"
        MS="Mean Square Error";
RUN;

PROC PRINT DATA=F NOOBS LABEL;
  WHERE Label1="Root MSE";
  VAR cValue1;
  LABEL cValue1="Root Mean Square Error";
RUN;

Meaningful Estimates (Linear Contrasts of Parameters)

The regression coefficients \(\beta_0, \beta_1, \beta_2\) in the quadratic regression model are difficult to interpret and not particularly meaningful in terms of understanding the relationship between mean serum total cholesterol and age. Quantities of more direct interest are mean serum cholesterol values for specific ages and differences in mean serum cholesterol values for two specific ages. For the quadratic regression model, the mean serum cholesterol value for 35 year olds is

\[\mu(35) = \beta_0 + \beta_1 \times 35 + \beta_2 \times 35^2 = \beta_0 + 35 \beta_1 + 1225 \beta_2 = \begin{pmatrix} 1 & 35 & 1225 \\ \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \end{pmatrix}\]

the mean serum cholesterol for 45 year olds is

\[\mu(45) = \beta_0 + \beta_1 \times 45 + \beta_2 \times 45^2 = \beta_0 + 45 \beta_1 + 2025 \beta_2 = \begin{pmatrix} 1 & 45 & 2025 \\ \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \end{pmatrix}\]

and the difference in mean serum cholesterol between 45 year olds and 35 year olds is

\[\mu(45) - \mu(35) = \beta_0 + 45 \beta_1 + 2025 \beta_2 - (\beta_0 + 35 \beta_1 + 1225 \beta_2) = 10 \beta_1 + 800 \beta_2 = \begin{pmatrix} 0 & 10 & 800 \\ \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \end{pmatrix}\]

Estimates of these quantities can be found by replacing the parameters \(\beta_0, \beta_1, \beta_2\) with the corresponding least squares estimates \(\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2\) in the above expressions. So, the estimated mean serum cholesterol value for 35 year olds is \[\hat{\mu}(35) = 3.1126 + 0.08122 (35) - 0.00075 (1225) = 5.0424 \] the estimated mean serum cholesterol for 45 year olds is \[\hat{\mu}(45) = 3.1126 + 0.08122 (45) - 0.00075 (2025) = 5.2585 \] and the estimated difference in mean serum cholesterol between 45 year olds and 35 year olds is \[\hat{\mu}(45) - \hat{\mu}(35) = 3.1126 (0) + 0.08122 (10) - 0.00075 (800) = 0.2160 \]

ODS EXCLUDE ALL;
PROC PLM RESTORE=QuadraticModel;
  ESTIMATE "Age 35" Intercept 1 AgePoly [1, 35];
  ESTIMATE "Age 45" Intercept 1 AgePoly [1, 45];
  ESTIMATE "Age 55" Intercept 1 AgePoly [1, 55];
  ESTIMATE "Age 65" Intercept 1 AgePoly [1, 65];
  ODS OUTPUT Estimates=E;
RUN;
ODS EXCLUDE NONE;

TITLE "Estimated Mean Total Cholesterol by Age";
PROC PRINT DATA=E NOOBS;
  VAR Label Estimate;
RUN;
TITLE;

ODS EXCLUDE ALL;
PROC PLM RESTORE=QuadraticModel;
  ESTIMATE "Age 45 vs. 35" AgePoly[1, 45] [-1, 35];
  ESTIMATE "Age 55 vs. 45" AgePoly[1, 55] [-1, 45];
  ESTIMATE "Age 65 vs. 55" AgePoly[1, 65] [-1, 55];
  ODS OUTPUT Estimates=E;
RUN;
ODS EXCLUDE NONE;

TITLE "Estimated Difference in Mean Total Cholesterol between Ages";
PROC PRINT DATA=E NOOBS;
  VAR Label Estimate;
RUN;
TITLE;

The estimated mean serum cholesterol levels are 5.04 mmol/L for 35 year olds, 5.26 mmol/L for 45 year olds (0.22 mmol/L greater than 35 year olds), 5.33 mmol/L for 55 year olds (0.07 mmol/L greater than 45 year olds) and 5.24 mmol/L for 65 year olds (0.08 mmol/L less than 55 year olds).

Centering and Scaling the Predictor Variable (Changing Units)

The intercept \(\beta_0\) in the quadratic regression model is the mean of the outcome with \(x=0\). If 0 is not a plausible value of \(x\) (or not in the range of the observed values of \(x\)), \(\beta_0\) does not represent a meaningful quantity and cannot be usefully interpreted. Interpretability can be improved by centering \(x\) to set the 0 point at a reasonable value within the range of \(x\). Often, the sample mean \(\bar{x}\) is used as a default centering point.

If we center \(x\) at \(x_0\), we can write the mean function of the quadratic regression model as \[\mu(x) = \gamma_0 + \gamma_1 (x-x_0) + \gamma_2 (x-x_0)^2\] instead of \(\mu(x) = \beta_0 + \beta_1 x + \beta_2 x^2\). After rearranging the terms in the centered model,

\[ \begin{eqnarray*} \mu(x) & = & \gamma_0 + \gamma_1 (x-x_0) + \gamma_2 (x-x_0)^2 \\ & = & \gamma_0 + \gamma_1 (x - x_0) + \gamma_2 (x^2 - 2 x x_0 + x_0^2) \\ & = & (\gamma_0 - \gamma_1 x_0 + \gamma_2 x_0^2) + (\gamma_1 - 2 x_0 \gamma_2) x + \gamma_2 x^2 \\ \end{eqnarray*} \]

it follows that \[\beta_0 = \gamma_0 - \gamma_1 x_0 + \gamma_2 x_0^2, \beta_1 = \gamma_1 - 2 x_0 \gamma_2, \beta_2 = \gamma_2\] or \[\gamma_0 = \beta_0 + \beta_1 x_0 + \beta_2 x_0^2, \gamma_1 = \beta_1 + 2 x_0 \beta_2, \gamma_2 = \beta_2\]

In our example, if we center the ages at 50 years, the estimated regression coefficients are \(\hat{\gamma}_0 = 5.3106, \hat{\gamma_1} = 0.006699, \hat{\gamma}_2 = -0.00075\) instead of \(\hat{\beta}_0 = 3.1126, \hat{\beta}_1 = 0.08122,\) and \(\hat{\beta}_2 = -0.00075\). Centering does not change the fitted values \(\hat{y}_i\), the residuals \(r_i\) or the residual variance \(\hat{\sigma}^2\). It is simply a tool for facilitating our interpretation of the model (and potentially making it easier to fit the model to our data).

DATA Analysis;
  SET Analysis;
    
  Age_C = Age-50;
RUN;

ODS EXCLUDE ALL;
PROC GLMSELECT DATA=Analysis;
  EFFECT AgePoly = POLYNOMIAL(Age_C / DEGREE=2);
  ODS OUTPUT ANOVA=A FitStatistics=F;
    
  MODEL TotChol = AgePoly / SELECTION=NONE;
    
  STORE QuadraticModel_C;
RUN;

PROC PLM RESTORE=QuadraticModel_C;
  ODS SELECT ParameterEstimates;
  SHOW Parameters;
RUN;
ODS EXCLUDE NONE;

DATA Score;
  DO Age = 20 TO 80 BY 1;
    Age_C = Age-50;
    OUTPUT;
  END;
RUN;

ODS EXCLUDE ALL;
PROC PLM RESTORE=QuadraticModel_C;
  SCORE DATA=Score OUT=Score Predicted LCLM UCLM LCL UCL / ALPHA=0.03125;
RUN;
ODS EXCLUDE NONE;

DATA All_C;
  SET Analysis Score(IN=S);
    
  Score = S;
RUN;

PROC SGPLOT DATA=All_C NOAUTOLEGEND;
  SCATTER X=Age Y=TotChol / MARKERATTRS=(SYMBOL=CircleFilled COLOR=LightBlue);
*  BAND X=Age Lower=LCL Upper=UCL / FILLATTRS=(COLOR=Red TRANSPARENCY=0.8);    
*  BAND X=Age Lower=LCLM Upper=UCLM / FILLATTRS=(COLOR=Red TRANSPARENCY=0.6);
  SERIES X=Age Y=Predicted / LINEATTRS=(COLOR=Red);
RUN;

PROC PRINT DATA=A NOOBS LABEL;
  WHERE Source="Error";
  VAR DF SS MS;
  LABEL DF="Residual df" 
        SS="Sum of Squares"
        MS="Mean Square Error";
RUN;

PROC PRINT DATA=F NOOBS LABEL;
  WHERE Label1="Root MSE";
  VAR cValue1;
  LABEL cValue1="Root Mean Square Error";
RUN;

ODS EXCLUDE ALL;
PROC PLM RESTORE=QuadraticModel_C;
  ESTIMATE "Age 35" Intercept 1 AgePoly [1, -15];
  ESTIMATE "Age 45" Intercept 1 AgePoly [1, -5];
  ESTIMATE "Age 55" Intercept 1 AgePoly [1, 5];
  ESTIMATE "Age 65" Intercept 1 AgePoly [1, 15];
  ODS OUTPUT Estimates=E;
RUN;
ODS EXCLUDE NONE;

TITLE "Estimated Mean Total Cholesterol by Age";
PROC PRINT DATA=E NOOBS;
  VAR Label Estimate;
RUN;
TITLE;

ODS EXCLUDE ALL;
PROC PLM RESTORE=QuadraticModel_C;
  ESTIMATE "Age 45 vs. 35" AgePoly[1, -5] [-1, -15];
  ESTIMATE "Age 55 vs. 45" AgePoly[1, 5] [-1, -5];
  ESTIMATE "Age 65 vs. 55" AgePoly[1, 15] [-1, 5];
  ODS OUTPUT Estimates=E;
RUN;
ODS EXCLUDE NONE;

TITLE "Estimated Difference in Mean Total Cholesterol between Ages";
PROC PRINT DATA=E NOOBS;
  VAR Label Estimate;
RUN;
TITLE;